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ABSTRACT 

We test four commonly used astrophysical simulation codes; Enzo, Flash, Gadget and 
Hydra, using a suite of numerical problems with analytic initial and final states. Sit- 
uations similar to the conditions of these tests, a Sod shock, a Sedov blast and both 
a static and translating King sphere occur commonly in astrophysics, where the ac- 
curate treatment of shocks, sound waves, supernovae explosions and collapsed haloes 
is a key condition for obtaining reliable validated simulations. We demonstrate that 
comparable results can be obtained for Lagrangian and Eulerian codes by requiring 
that approximately one particle exists per grid cell in the region of interest. We con- 
clude that adaptive Eulerian codes, with their ability to place refinements in regions 
of rapidly changing density, are well suited to problems where physical processes are 
related to such changes. Lagrangian methods, on the other hand, are well suited to 
problems where large density contrasts occur and the physics is related to the local 
density itself rather than the local density gradient. 
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1 INTRODUCTION 

' Over the last decade, numerical simulations have devel- 
oped into a cornerstone of astrophysical research. From in- 
teractions of vast clusters of galaxies to the formation of 
proto-planetary discs, simulations allow us to evolve sys- 
tems through time and view them at every angle. They are 
an irreplaceable test-bed of our physical understanding of 
the Universe. 

A number of hydrodynamics codes have been developed 
that are widely used in this field and still more are being de- 
veloped. Fundamentally, they all do the same job; they solve 
the equations of motion to calculate the evolution of matter 
through time. Whether the matter represents a nebula for 
the birth of a star or a network of galaxy clusters, the basic 
technique remains the same. 

However, the algorithms used to solve these equations 
vary from code to code and this results in differences in the 
resulting data. Understanding the origin of these variations 
is vital to the understanding of the results themselves; is an 
observed anomaly an interesting piece of new physics or a 
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numerical effect? As observational data takes us deeper into 
the Universe, it becomes more important to pin down the 
origin of these numerical artifacts. 

Additionally, it is difficult to compare results from sim- 
ulations run with different codes. With observations, papers 
clearly state the properties of the instrument such as the 
diameter of the mirror and the wavelengths it is most sen- 
sitive to. While a brief description of the code is always in- 
cluded in theoretical papers, there exists no obvious conver- 
sion to other numerical techniques and therefore the results 
are more difficult for the reader to interpret. 

The problem of code comparison is not new and it is a 
topic that has recently created a great deal of interest. The 
reason for its current importance is a positive one; improved 
numerical techniques and increased computer power have re- 
sulted in simulations reaching greater resolutions than could 
have been imagined even a few years ago. However, this high 
refinement comes at a price; as we start to pick out the de- 
tail of these complex fluid flows, the physics we need to con- 
sider gets dramatically more complicated. This brings us to 
the main question code comparison projects are trying to 
answer; can we use the same tools for this new regime of 
problems? 
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A number of papers have come out that tackle this. One 
of the most famous is the Santa Barbara Comparison Project 
ijFrenk et al.l 1 19991 ) which compared many of the progeni- 
tors of today's codes by running a model of a g a laxy c luster 
forming. Taking a different tack iThacker et al] (I2OO0I ) com- 
pared the performance of a dozen different implementations 
of a single approach (in this case smoothed particle hydro- 
dynamics, SPH) on standard astrophysical problems that 
included the Sod shock, examining the range of outcomes 
that were available from a single technique. They concluded 
that one of the weaknesses of SPH was the weak theoreti- 
cal grounding which allows several equally vi able formula- 
tions to be derived. More recent work includes I Agertz et al.l 
who focus specifically on the formation of fluid insta- 
bilities, comparing the formation of Kelvin-Helmholtz and 
Rayleigh- Taylor instabilities in six of the most utilised codes. 
Additionally, O'Shea ct al. (2005) has completed a direct 
comparison between two particular codes [Enzo and Gad- 
get2, see below) looking at the formation of galaxies in a 
cosmological context. 

All these projects give detailed insights into the differ- 
ences between the codes, but are unable to provide a quan- 
titative measure of how well a code performs in a particu- 
lar aspe ct. This is e s pecial l y true of the c o smolo gical-based 
tests of iFrenk et al] (|l999D : lO'Shea et all Hooi) where the 
problem is not sufficiently well-pos ed for convergence onto a 
single answer. lAgertz et al.l (|2007t ) sets out a simpler prob- 
lem and compares it to analytical predictions, but the sys- 
tem is still sufficiently complex not to have an exact solu- 
tion. Additionally, no previous comparison has attempted to 
quantitatively compare different codes to one another; ask- 
ing whether it is possible to obtain identical results and with 
what conditions. Without this crucial piece of information, 
it is impossible to fully assess a piece of work performed 
by an unfamiliar code or to judge which code might be the 
most suited to a given problem type. This has resulted in 
somewhat general comments being made about the differ- 
ences between numerical techniques which has led to many 
myths about a code's ability becoming accepted dogma. 

The set of tests we present in this paper are designed to 
tackle these difficulties with the intention that they might 
become part of an established test programme all hydro- 
dynamical codes should attempt. We present four problems 
that specifically address different aspects of the numerical 
code all of which have expected 'correct' answers to compare 
to. The first two tests, the Sod shock tube and Sedov blast, 
are both strong shock tests with analytical solutions. The 
third and forth tests concern the stability of a galaxy clus- 
ter and are primarily tests of the code's gravitational solver. 
For all four tests we directly compare the codes against the 
analytic solution and present an estimate of the main sources 
of any systematic error. 

The remainder of this paper is organised as follows: in 
section 2 we give a short summary of the main features of 
each of the four codes we have employed. In section 3 we 
deal with the Sod shock and Sedov blast tests, setting out 
the initial and final states and comparing each code against 
them. We repeat this exercise for a static and translating 
King sphere in section 4. Finally we discuss our results and 
summarise our conclusions in section 5. 



2 DESCRIPTION OF THE CODES 

The two major techniques for modelling gases in astro- 
physics are smoothed particle hydrodynamics (SPH) and 
adaptive mesh refinement (AMR). In the first of these, the 
gas is treated as a series of particles whose motion is dictated 
by Lagrangian dynamics. In AMR, the gas is modelled by 
a series of hierarchical meshes and the flow of material be- 
tween cells is calculated to determine its evolution. There 
are a variety of codes which utilise both techniques and four 
of the major ones will be used to run the tests presented in 
the this paper, two of which use SPH {Hydra and Gadget2) 
and two of which use AMR {Enzo and Flash). 

2.1 Enzo 

Enzo is a massively parallel, Eulerian adaptive mesh re- 
finemen t code, capable of both hydro and N-bod y calcu- 
lations (|0'Shea et al.l 12004 iBrvan fc Normanlll997l ). It has 
two hydro- algorithms which can be selected by the user; the 
piecewise parabolic method (PPM) and the Zeus astrophys- 
ical code. The PPM solver (Woodward & ColcUa 1984) uses 
Godunov's method but with a higher-order spatial inter- 
polation, making it third-order accurate. It is particularly 
good at shock capturing and outfiows. The Zeus method in 
Enzo is a three-dimensional implemen tation of the Zeu s as- 
trophysical code developed bv .Stone fc NormanI (|l992l fl. It 
is a simple, fast algorithm that allows large problems to be 
run at high resolution. Rather than Godunov's method, Zeus 
uses an artificial viscosity term to model shocks, a technique 
which inevitably causes some dissipation of the shock front. 
We compare both these hydro-schemes in these tests. 

2.2 Gadget2 

Gadgets is a massively parallel, Lagrangian, cosmological 
code that is publicly available from the author's website 
(Springcl 2005J. It is an N-body/SPH code that calcu- 
lates gravitational f orces by means of the Tree method 
(|Barnes fc Hut|[l986l ) and is also able to optionally employ 
a Tree-PM scheme to calculate the long range component 
of the gravitational interactions. In order to follow the hy- 
drodynamic behaviour of a coUisional medium, the code 
us es the entropy-conserving fo rmulation of SPH described 
in ISpringel fc HernquistI (|2003l ) : the maii i difference of this 
approach with respect to the standard iMonaghanI (|l992l ) 
formulation of SPH resides in the choice of describing the 
thermodynamic state of a fluid element in terms of its spe- 
cific entropy rather than its specific thermal energy. This 
leads to a tight conservation of both energy and entropy in 
simulating dissipation-free systems. Additionally, Gadget2 
employs a slightly modified parametrisation of the artificial 
viscosity (by intr oducing the so called "signal velocity" as in 
lMonaghanl ( |200ll )'). The user is allowed to set the strength of 
this artificial viscosity for the specific problem being consid- 
ered via an input parameter, ArtBulkVisc. The time step- 
ping scheme adopted by the code is a leap-frog integrator 

^ Note that the hydrodynamics code ZEUS-3D has been devel- 
oped independently of Enzo (Zeus) and its performace is not 
equivalent. 
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guaranteed to be symplectic if a constant timestep for all 
particles is employed. In this work we have exploited the 
possibility of using fully adaptive individual timesteps for 
all the particles in the simulation, this being a standard 
practice. 



2.3 Hydra 

Hydra is an adaptive particle-particle, particle-mesh 
code combined with smoothed pa rticle hydrodynamics 
llCouchman. Thomas, fc Pearcd liggsll. It has the significant 
disad vantage that even though m assively parallel versions 
exist (jPearce fc CouchmanI 119971 ) , the publically available 
version is not a parallel implementation and so this code 
cannot, as released, be used for very large simulations. Akin 
to Gadgets, Hydra uses an entropy conserving implemen- 
tation of SPH, but unlike Gadget2, Hydra does not have 
fully adaptive individual timesteps. Although the timestep 
adapts automatically from one step to the next all the par- 
ticles move in lockstep. 



2.4 Flash 

Flash is a publicly available massively parallel Eulerian 
AMR code developed by the A lliances Center for Astro- 
physical Thermonuclear Flashes (|Frvxell et al.llioool '). Origi- 
nally intended for the study of X-ray bursts and supernovae, 
it has since been adapted for many astrophysical condi- 
tions and now includes modules for relativistic hydrodynam- 
ics, thermal conduction, radiative cooling, magnetohydro- 
dynamics, thermonuclear burning, self-gravity and particle 
dynamics via a particle-mesh approach. Flash uses the oct- 
tree refinement scheme of the PARAMESH package, with 
each mesh block containing the same number of internal 
zones. Neighbouring blocks may only differ by one level of 
refinement with each level of refinement changing the reso- 
lution by a factor of two. The hydrodynamics are based on 
the PROMETHEUS code (FryxcU, Miillor, fc Arn ctt 1989|). 
The input states for the Rieman n solver are obtained using a 
directionally split PPM solver (|Woodward fc Colellalll984l ') 
and a variable time step leapfrog inte grator with second 
order Strang time splitting is adopted jstrang| 'l968'). This 
work uses a modified hybrid FFTW based multigrid solver 
to solve Poisson's equation and determine the gravitational 
potential at each timestep. This results in a vast reduction 
in time spent calculating the self-gravity of the simulation 
relative to a conventional multigrid solver. Flash's refine- 
ment and de-r efinement criteria can incorporate the adapted 
iLohnej (|l987l ) error estimator. This calculates the modified 
second derivative of the desired variable, normalised by the 
average of its gradient over one cell. 



3 SHOCKING TESTS 

One of the greatest differences between simulations per- 
formed now versus those undertaken five years ago is the in- 
creasing importance of modelling strong shocks accurately. 
While it has long been known that the Universe is a vi- 
olent place, with events such as supernovae, galaxy merg- 
ers and AGN generating blasts which rip through the in- 
tergalactic medium, simulations did not have the resolution 



Figure 1. Example projections of the initial conditions for a 
three-dimensional Sod shock test. Black and white regions repre- 
sent fluids of different densities. The left-hand image shows the 
shock face oriented along the [1,0,0] plane, while the right-hand 
image shows it oriented along the [1,1,0] plane. In actuality our 
second test is oriented in the [1,1,1] plane i.e. oblique to all the 
axes. 



to see such phenomena in detail, so these sharp disconti- 
nuities were largely ignored. Now, as we struggle to un- 
derstand the effects of feedback in galaxy formation , mul- 



tipha se media are essential physics ( Tasker fc Brva nll2008l , 
20061; IWada fc Normanll2007l : iRobertson fc Kravtsovlf2007l ). 

In order to attack such problems codes must be able to cap- 
ture shocks with some proficiency. These two problems, the 
Sod shock test and the Sedov blast test, explicitly test the 
resolution of shock jumps and allow comparison with exact 
analytical solutions. 



3.1 Riemann Shock Tube Problem 



The shock tube problem l|Sodlll978l ) has been used exten- 
sively to test the abilit y of hydrodynam i cs codes to resolve a 
sharp shock inte rface (|Feng et al.ll200 IShapiro et al. I ll996l : 
IRvu et al.lll993l ). The test set-up is simple, consisting of 
two fluids of different densities and pressures separated by a 
membrane that is then removed. The resulting solution has 
the advantage of showing all three types of fluid discontinu- 
ities; a shock wave moving from the high density fluid to the 
low density one, a rarefaction (sound) wave moving in the 
opposite direction and a contact discontinuity which marks 
the current location of the interface between the two fluids. 
For this test the initial conditions are traditionally chosen 
such that the pressure does not jump across the contact dis- 
continuity. 

We extend the traditional one-dimensional shock tube 
problem to consider two three-dimensional set-ups; the first 
of these has the fluid membrane at 90° to the x-axis of the 
box ([1,0,0] plane), causing the shock to propagate parallel 
to this axis. In the second test, the membrane is lined up 
at 45° to each of the x,y and z axes ([1,1,1] plane). This 
change in orientation of the shock is designed to highlight 
any directional dependencies inherent in the code, as illus- 
trated in Figure [1] All analysis was performed perpendicular 
to the original shock plane. 

For our particular set-up we chose the initial density 
and pressure jump either side of the membrane to be from 
(pi = 4,pi = 1) to (p2 = 1,P2 = 0.1795) with the fluid 
initially at rest. The poly tropic index was 7 = |. Periodic 
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boundary conditions were used and the results were analysed 
at t = 0.12. 



3.1.1 General Results 

Figures [2] and |3] show the results from all four codes running 
this test. In the results presented in these figures, both Enzo 
and Flash used an 100'' initial (minimum refinement) gricjf] 
with two levels of higher refinement each of which decreased 
the cell size by a factor of 2. The smallest cell size in this 
case was therefore 0.0025 of the unit box. Enzo refined any- 
where the gradient of the derived quantities exceeded a crit- 
ical value whe reas Flash placed refinements according to the 
lLohnej (j 19871) error estimator described in section [2^ which 
places subgrids based on the second derivative of the derived 
quantities. For the SPH codes, both Hydra and Gadgets were 
run with 1 million particles formed from two glasses contain- 
ing 1.6 million and 400,000 particles. The solid black line in 
all cases is the analytical solution which requires a shock at 
X = —0.095 and a contact discontinuity in the density at 
X = -0.032 at t = 0.12. 

What is clear from Figures [2] and [3] is that all of the 
codes pass the zeroth level test and successfully reproduce 
the shock jump conditions, although both the SPH codes 
suffer from visible ringing and broadening around any dis- 
continuities. Enzo (Zeus) does not produce the shock jump 
condition as accurately as Enzo (PPM) and Flash, as seen in 
the plots of energy and entropy in Figure [S] where its post- 
shock values are around 2% lower. In the oblique case, the 
quadratic viscosity term, QV, in Enzo (Zeus) was increased 
from its default value of 2.0 to 10.0. The effect of this value 
is discussed more fully in relation to the Sedov blast test 
in section 13.2.31 For the planar case, QV was kept at 2.0. 
Pleasingly none of the other codes appear to have any visi- 
ble directional dependence, performing equally well in both 
the grid aligned and oblique cases we tried. All the codes 
could equally well resolve the location and smooth rise of 
the rarefaction wave but both the SPH codes struggle with 
the contact discontinuity, with a large overshoot not seen by 
either of the mesh based codes. This is partly due to the ini- 
tial conditions as for the SPH codes the sudden appearance 
of a density jump introduces a local source of entropy. In 
this paper we are contrasting the results from the different 
approaches rather than studying the Sod shock problem it- 
self in detail. As it is a standard test case higher resolution 
resu lts can be found in the individual code's method paper s 
(e.g. iFrvxell et alll2000l : iThacker et alll2000l : ISpringelllioOsI ) . 

Closer inspection of the data reveals differences in each 
code's capacity to handle strong shocks. In the bottom row 
of Figure [21 we show a close-up of the density over the shock- 
front. Capturing shocks accurately is an area that SPH codes 
traditionally struggle with more than their Eulerian coun- 
terparts due to their inherent nature of smoothing between 
particles. Indeed, we see in this figure both Gadget2 and 
Hydra have a smeared out the interface compared to Enzo 
and Flash 's steep drop in density. Small differences between 

^ In Enzo's case, this is the 'root' grid. Flash does not have a root 
grid, but can specify a minimum cell size that must be maintained 
everywhere. For ease of typing, we shall refer to the coarsest grid 
in both cases as the 'initial' grid. 



the AMR codes are also visible here. Enzo (PPM) spreads 
the shock front over three cells, whereas Enzo (Zeus's) use 
of the artificial viscosity term extends this to five. Flash's 
PPM scheme gives very similar results to Enzo (PPM), also 
spreading the shock front over three cells. 

Figure [3] shows the pressure, internal energy, velocity 
and computational entropy = — ^rr^ over the region of 
interest for the planar [1,0,0] set-up (left column) and the 
oblique [1,1,1] set-up (right column). 

Both Hydra and Gadgets show signs of post-shock ring- 
ing in the velocity plot. The two lagrangian codes adopt 
different implementations of artificial viscosity; with a fre- 
quently used choice of the viscosity parameter (ArtBulkVisc 
= 1) the ringing features in Gadget2's profiles appears to be 
more pronounced than in Hydra's (not shown in the plot). 
In order for Gadget2 to get closer to Hydra's performance, a 
choice of a significantly higher viscosity parameter has been 
necessary (namely ArtBulkVisc = 2). The results produced 
under such a choice are shown in Figure [3] The SPH codes 
also exhibit a large spike in both internal energy and en- 
tropy at the location of the contact discontinuity. This is 
driven by the initial conditions, where two independent par- 
ticle distributions suddenly appear immediately adjacent to 
one another. 



3.1.2 Quantitative Comparison 

We saw in the previous section that most of the codes model 
shock development and sound wave propagation with rea- 
sonable success and that they show no directional preference 
to the orientation of the shock interface. However, differences 
were apparent between each of the codes, most obviously 
between the SPH and AMR techniques (unsurprising, since 
their numerical algorithms fundamentally differ). In this sec- 
tion, we quantitatively compare results from each code and 
attempt to get as close a match between their results as 
possible. For simplicity, we confine the AMR codes to using 
static meshes for this comparison. 

Figure |4] shows a graphical comparison of the density 
projection between a Hydra run of 1 million particles and 
Enzo (PPM) for different grid sizes. It is clear that a grid size 
of 20^ produces significantly poorer results than the Hydra 
data whereas a grid size of 100^ produces significantly bet- 
ter results particularly in the low density region. Although 
the situation is confused by the lack of points available with 
Enzo, a grid size of 20^ can produce no more than 20 distinct 
values across the length of the volume being modelled, it is 
clear that this is too few to recover this model. However, even 
with 50^ cells in the box the major features are largely recov- 
ered. To make further progress we require a more quantita- 
tive way of comparing the results from the codes. To achieve 
this, we employ a cubic spline to interpolate the data from 
all runs to the same 178 x points at which we have calcu- 
lated the analytical values. The residue between each new 
curve and the analytical solution was then summed and di- 
vided by the number of points. We considered the residue 
both across the whole region of interest from [-0.15, 0.15] 
and just across the shock-front from [-0.13, -0.06]. 

The top part of table [T] shows the residues from Enzo 
over these ranges for increasing static grid size. As Figure |4] 
showed, the overall fit across the whole region of interest 
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Figure 2. Projected density parallel to the shock front from the 3-dimensional shock tube test. Left-hand plots show results from the 
[1,0,0] shock set up, right-hand from the [1,1,1] set-up. Both are at t = 0.12. The bottom plots show a close-up of the shock front itself. 



improves rapidly with the number of grid points. The im- 
provement for the AMR run highlights one of the greatest 
assets of adaptive grid coding; the ability to tag cells for re- 
finement based on slopes or shock-fronts, allowing the extra 
resolution to be concentrated on these problematic areas. 
Table [1] also shows that Enzo ( Zeus ) typically has poorer 
results than Enzo (PPM) for equal numbers of grid cells, 
something that came out in Figures [5] and |31 and is a result 
of energy loss due to the inaccurate treatment of the shock 
jump conditions by Enzo (Zeus). 

The lower part of table [1] shows the residues for Hydra 



and Gadget2 at different resolutions. Both SPH codes show 
similar results, with similar residues in each at the same res- 
olution. If we compare these values to the nearest Enzo result 
we can get an estimate for the number of particles needed 
per grid cell in the low density region for similar recovery of 
the shock profiles. Included after the number of particles for 
each of these runs is the effective number of particles in the 
high and low density regions. As is obvious from this table 
and Figure |4] roughly matching accuracy occurs somewhere 
between these two limits. This is equivalent to when there 
is approximately one SPH particle per grid cell. 




Figure 3. Projected data parallel to the shock front from the 3-dimensional shock tube test. The left hand column shows the results from 
the [1,0,0] shock front alignment, while the right hand column gives those from the [1,1,1] set-up. Top-to-bottom arc plotted pressure, 
internal energy, velocity and entropy. As with figure|2]all the panels are from t = 0.12. 
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Figure 4. Comparison between data from Enzo (PPM) at different static grid sizes with Hydra's run of 1 million particles. Results 
plotted show the projection of the density for the [1,0,0] set-up. 



Table 1 . Residue of from the analytical solution of the Sod shock 
density in the planar [1,0,0] set-up for different static grid sizes 
and SPH resolutions for Enzo, Hydra and Gadget2. 



Enzo (PPM) residue Enzo (Zeus) residue 
Grid size Interface (Shock) Interface (Shock) 



10^ 0.343 (0.383) 0.238 (0.220) 

20^ 0.232 (0.166) 0.173 (0.198) 

303 0.164 (0.164) 0.165 (0.172) 

40^ 0.129 (0.125) 0.166 (0.174) 

50^ 0.095 (0.089) 0.119 (0.130) 

60^ 0.094 (0.110) 0.114 (0.118) 

70^ 0.076 (0.082) 0.105 (0.126) 

80^ 0.086 (0.082) 0.081 (0.090) 

90^ 0.061 (0.075) 0.088 (0.088) 

100^ 0.069 (0.064) 0.085 (0.079) 

AMR run (400^) 0.032 (0.030) 0.035 (0.036)) 



Hydra residue Gadget2 residue 
No. particles Interface (Shock) Interface (Shock) 



8M (2343,148^) 0.0559 (0.065) 
IM (117^,743) 0.0852 (0.104) 0.0817 (0.089) 

250k (743, 48^) 0.112(0.121) 0.113(0.102) 



The analysis done in this section does have one obvi- 
ous flaw; introducing adaptive meshes allows the grid-based 
codes to achieve a high resolution while employing many less 
cells overall. Our equivalence of one SPH particle per cell 
applies to cells in the high resolution region not the total 
number of cells in the simulation. As such, this comparison 
works best in situations where the resolution is needed in 
the highest density areas. In both this problem type and the 
one described in the next section, the high density areas are 
not those that require higher resolution, rather the region of 
interest is that where the density is changing rapidly. Since 
usually SPH particles are tied to the mass flow of a system. 



they are not able to increase the resolution over shock fronts, 
unlike their grid counterparts. However, in many astrophys- 
ical problems, such as those we shall meet in section U] this 
is not necessarily a serious limitation. 

3.2 The Sedov Blast Wave Test 

The Sedov Blast Test (|Sedov|[l959l ) is an intense explosion 
caused by a quantity of energy deposited in the centre of the 
simulation box. The result is a strong spherical shock that 
propagates through the background homogenous medium. 

This test is particularly appropriate in astrophysics 
since it represents well the physics required to deal with 
supernova explosions. An inability to resolve the resulting 
shock front will lead to an incorrect estimate of the energy 
deposited into the galaxy and the volume it affects. 

It also poses different problems to both the particle and 
grid schemes described above. Although the shock front ex- 
pands as it travels outwards it also sweeps up particles which 
increase the density contrast of the shell relative to the am- 
bient medium. The particle codes will therefore find it pro- 
gressively easier to reach the required resolution. Since the 
shock front is also spherical, grid codes will have to cope 
with any artiflcial grid-alignment effects. It is therefore a 
test that is both appropriate and taxing. 

The analytical c alculat i on of the sh ock propagation is 
described in full in ISedovl { I959I ) and iLandau fc Lifshit3 
(jl959i ) who show the shock front's radius is given by: 

Kt)=f-^)''%^/^ (1) 

where Eo is the initial energy injected, po is the background 
density and a = 0.49 for an ideal gas with 7 = 5/3. 

Since there is no cooling or gravity in this test, the 
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Figure 5. Density projections of the Sedov blast test performed using Enzo (PPM) at t = 0.025, 0.05, 0.075 and 0.1. 



quantities above are all unitless. In our calculations, the ini- 
tial density of the medium is po = 1, the explosion energy 
Eo = 10^. The shock expands into a box of side 10 and 
we compare our results at t — 0.1. Figure [5] shows density 
projections (density integrated along the length of the sim- 
ulation box) of the simulation box during the Sedov test for 
times from 0.025 to 0.1. 

3.2.1 General Results 

Figure |6] shows the results from the Sedov blast test at 
t = 0.1 for each code and the analytical solution. For this 
test, Enzo and Flash used an initial grid of 100'^ with two 
levels of refinement, each reducing the cell size by a factor of 
two. These subgrids were placed where the shock front was 
detected. This gave a minimum cell size of 0.025 over the 
shock. For the AMR codes the central energy was injected 
into a spherical region of radius 3.5 grid cells on the finest 
levels, corresponding to r = 0.0875, a common choice for 
this test. Both Hydra and Gadgets used 1 million particles 
relaxed into a glass, with the energy added to the central 32 
particles in a top hat distribution. 

We have deliberately chosen a large energy jump for 
our Sedov blast in order to make this test challenging. A 
short discussion is added to the end of this section describ- 
ing some of the pit-falls we encountered while attempting to 
successfully complete it. What is immediately apparent from 
Figure [H] is that Enzo (Zeus) performs significantly worse on 
this test than all the other codes, with its shock edge lag- 
ging substantially behind the analytical position. It shows 
a significant energy loss in the first few timesteps, loosing 
25% of its internal energy during this stage. In the previ- 
ous section, we saw that Enzo (Zeus) spread the shock front 
out slightly more than Enzo (PPM) and Flash, which both 
utilise Godunov schemes, but the difference was minimal and 
the shock-front was still sharp and well represented. Here, 
however, we see that Enzo (Zeus) underestimates the shock 
position at t = 0.1 by almost 4%. The cause of this sig- 
nificant energy loss is the production of a diamond-shaped, 
rather than spherical, shock front at very early times. This 
is discussed further in section [3. 2. 3 1 Plots of the internal en- 
ergy and pressure in Figure[S]also add to illustrate this effect, 
showing the eff'ect of Enzo (Zeus)^s early energy loss and the 
resulting low pressure prior to the shock front. Both Enzo 
(PPM) and Flash perform this test extremely well, match- 
ing the analytical solution to the shock's edge. The two SPH 
codes successfully recover the location of the peak of the Se- 



dov blast shell but it is smoothed out in the radial direction, 
producing a lower peak density and a broader shell. For the 
Gadgets run the gas velocity interior to the blast appears 
to be enhanced in Figure [6] This is an artifact of the en- 
tropy scatter discussed below. The gas interior to the blast 
front will rapidly sort itself in entropy, resulting in smooth 
density and temperature profiles, leaving the velocity profile 
disordered (and so enhanced). 

Figure [7] shows density projections of the final state of 
the simulations at t — 0.1. The analytical position of the 
shock front is shown by the solid red line. As can be seen 
on Figure [7] although the location of the shock front in the 
Gadgets figure is well recovered there seems to be a lot of 
"noise" producing a grainy appearance to the shaded image. 
Although energy and entropy are well conserved by Gadgets 
for extreme shocks like this one the Gadgets viscosity im- 
plementation introduces a lot of entropy scatter. The high 
entropy particles are too hot for their surroundings and drive 
small bubbles, producing the structure seen in Figure [T] This 
feature of Gadgets does not appear for more normal shock 
jumps such as that studied in the previous section or for less 
energetic Sedov blasts. We have tried several approaches to 
resolve this problem which are mentioned in section 13.2.31 
below for those interested. 



3.2.S Quantitative Comparison 

As the previous section showed, four of our codes appear 
to do a reasonable job of following the location of the the 
Sedov blast front. Due to its inability to conserve energy 
initially, Enzo (Zeus) does not perform this test well and 
only achieves its current results by adapting code parame- 
ters (see "code-specific issues" below). In order to quantita- 
tively compare the codes we have measured the position of 
the peak of the density, the maximum density obtained and 
the width of the shock front at half this maximum density 
as a function of time. The results of this exercise are given in 
Figure ISl We see that Enzo (PPM), Flash, Hydra and Gad- 
gets are all capable of accurately following the location of 
the shock front. As is visible on Figures[6]and[8]the two SPH 
codes do not obtain the high maximum densities recovered 
by the AMR codes. GadgetS's maximum peak position is 
lower than Enzo (PPM) and Flash by 8% and Hydra is 13% 
less. Enzo (Zeus)'s lag in the shock front propogation makes 
its position equivalent to a time of t = 0.085, at which time 
it is between the PPM-based schemes and the SPH tech- 
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Figure 6. Results from the Sodov blast test for each code at t = 0.1. Clockwise from top left shows density, velocity, pressure and 
internal energy. The black solid line marks the analytical solution. 




Figure 7. Density projections of the Sedov blast test at t = 0.1 for all codes, with the expected shock position marked by the red circle. 
Left to right shows Enzo (PPM), Enzo (Zeus), Flash, Gadget2 and Hydra. The density ranges from [10°-^, 101-23]. 



this disadvantage, the SPH codes correctly reproduce the 
shock front position and, at around 10%, the differences in 
their peak density values are remarkably small. 

It is difficult to undertake a "number of particles per 
grid cell" comparison similar to the one above for this test 
for two reasons: firstly, the initial conditions rapidly become 
impossible to set up in a spherically symmetric way for the 
grid codes if they are not allowed to refine heavily. Our 
initial condition set-up inserts the energy into a region of 



niques with a top density value 6% below Enzo (PPM) and 
Flash. 

The lower SPH peak density is not a surprise as the in- 
herent smoothing of the SPH method broadens the shock 
front and lowers the maximum recoverable density. The 
AMR codes have the ability to insert additional resolution 
elements where quantities are rapidly changing rather than 
where lots of mass has piled up and this acts to improve their 
fit to the analytic solution in the region of the blast. Despite 
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Figure 8. For the Sedov blast test, from left to right; position of shock front over time, the maximum density value and the width of 
the shock front at half maximum. 



radius 3.5 grid cells (or 180 cells). If fewer than this num- 
ber are used the cells do not represent a spherical energy 
input well. Without refinement a region encompassing this 
many cells rapidly becomes large, whereupon the initial con- 
ditions no longer represent a Sedov explosion. In addition, 
small numbers of large cells simply cannot be expected to 
follow a spherical blast and little useful is learnt. Secondly, 
the SPH codes are helped at late times because the Sedov 
blast "sweeps up" the particles it encounters (as it is sup- 
posed to) which effectively increases the number of SPH 
resolution elements in the shock front, as is to be expected 
given the density enhancement. The AMR codes naturally 
account for this by adding extra layers of refinement. How- 
ever, if an unrefined comparison is attempted this high level 
of refinement has to be present everywhere and at all times, 
so an incredible number of cells are needed. For this partic- 
ular test we find that we need static grids of order 250'' to 
reproduce a blast with similar resolution to an SPH model 
with only lOO"' particles. The quantity that should be di- 
rectly compared is the number of cells per particle in the 
region of interest. As above we find that roughly one SPH 
particle per AMR cell is required to obtain similar effective 
resolutions. 



3.2.3 Code-Specific Issues for the Sedov blast 

Figure |9]shows the density projection of the Sedov blast test 
performed by Enzo (Zeus) at early times. The problem be- 
comes immediately apparent; instead of being spherical, the 
shock-front is an assymetrical diamond shape. Over time 
the shock-front becomes spherical, as can be seen in Fig- 
ure [T] and the dramatic energy loss stops, but by this stage 
25% of the shock's energy has been lost. The reason Enzo 
(Zeus) fails to produce a spherical shock-front at early times 
is not obvious. It appears overly sensitive to the initial set- 
up not being a perfect Sedov blast start (since a point-like 
energy injection can only be approximate). Many difi'erent 
solutions were attempted, including injecting the energy as 
a Gaussian rather than 'Top hat' profile and adding it into 
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Figure 9. Density projection showing the early evolution of the 
Sedov blast test at t = 0.01 for the AMR codes. Left to right 
shows Enzo (PPM), Enzo (Zeus) and Flash. The projected den- 
sity range is [lO'^'*'^, 10^'^]. At this time, the assymmetry of the 
shock- front with Enzo (Zeus) is clearly visible. Both SPH codes 
give spherically symmetric results. 



a larger radius (which helped, but did not greatly improve 
the situation). Even to recover the pretty miserable results 
that Enzo (Zeus) achieved above we had to use a somewhat 
dirty fix of increasing the code's quadratic artificial viscosity 
term, QV. This value is a broadening parameter that con- 
trols how many cells the shock is spread over. By default 
QV — 2, wh i ch is used for the tests in sections 13.11 and U 
lAgertz et all (|2007l ) found that varying QV made little dif- 
ference to the evolution of the fluid, but just broadened the 
shock. In this case, we find that increasing the width of the 
shock decreases the assymmetry at early times, although 
Figure shows the effect was far from perfect. Using the 
default value of QV = 2, the shock-front lags even further 
behind the result shown in Figure [S] which uses QV = 10. 
The increase in QV comes at a price; the spreading of the 
shock over more cells weakens it, causing the peak density to 
be lowered. Raising QV beyond 10 improves the position of 
the shock further, but the value now becomes unphysically 
high and we do not recommend using it. 

The Sedov bl as t ha s already been used by 
ISpringel fc Hernguistl (|2002l ) to test Gadqet2, an d in - 
deed for the parameters ISpringel fc Hernguistl (|2002l ) 
used Gadget2 works beautifully and does not exhibit any 
spurious entropy driven bubbles. In order to investigate 
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their origin we tried injecting the energy into a single 
point particle or smoothing the injection profile using a 
Gaussian function (with a cut-off at a few percent of its 
central value) rather than a top hat, and also changing 
the number of SPH neighbours in order to help numerical 
convergency. None of these changes makes any practical 
difference to our test and the entropy driven bubbles 
still occurred. In all the Gadget2 runs performed for this 
test, the artificial viscosity parameter ArtBulkVisc has 
been set to the frequently used value of 1; there are hints 
that a choice for a higher value might lead to a sensible 
reduction of the entropy driven bubbles (see Pakmor et 
al., in preparation). It appears that in extreme shocks the 
viscosity implementation employed by Gadget2 does not 
sufficiently prevent particle inter-penetration and we warn 
users of Gadgets to be wary of the generation of spurious 
entropy scatter in the vicinity of extreme shocks. 



4 GRAVITATIONAL TESTS 

In this section, we move away from the formation and resolu- 
tion of shocks to look at a new aspect of the codes; how they 
deal with gravity. While treatment of fiuids is important, few 
astrophysical simulations can be performed without a self- 
gravitating gas. However, adding self-gravity, where every 
fluid element is affected by every other, dramatically com- 
plicates the situation and it is not possible to design a test 
with an exact analytical solution anymore. Since it is still 
essential for the purpose of this comparison that our prob- 
lems remain well-posed, we select situations in which the 
correct behaviour of the system is known, even if it cannot 
be mathematically expressed. To do this, we perform two 
tests; the first of these concerns a static gas profile in equilib- 
rium. Gravity acts to try and collapse the gas, while pressure 
pushes it outwards. While these forces remain perfectly bal- 
anced, the gas remains at rest. This situation is analogous to 
a relaxed galaxy cluster and requires the code to resolve the 
gas density over many orders of magnitude. The second test 
involves the same cluster translating through the box. By 
using periodic boundary conditions, the cluster's velocity is 
chosen such that it should return to its original position af- 
ter 1 Gyr. With no external forces, the cluster should remain 
in hydrostatic equilibrium and retained its profile during the 
translation. 



4.1 Initial conditions for the cluster 

The model used for the galaxy cluster is the King model 
ifKhi g 1966; Padmanabhan 2002), which was chosen because 
it possesses a finite radial cut-off, and is therefore a well 
defined problem for a code comparison. Its form is based on 
the distribution function: 
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where e = — , is the coordinate change for the shifted 
energy, pc is the central density and a is related to (but 
not equal to) the velocity dispersion. The resulting density 
distribution of this cluster vanishes at the tidal radius, rt . 
Integrating over all velocities yields a density distribution: 
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Putting this into the Poisson equation results in a second or- 
der ODE which can be solved numerically. This model has 
three independent parameters, the mass of the cluster, the 
tidal radius and the concentration c = log^Q(rt/ro), where 
ro is the central or King radius. For this problem, we se- 
lected a concentration of 3, vt — 1 Mpc and a cluster mass 
of 10^'' M0. This results in a King radius of 1 kpc. Therefore, 
to successfully maintain hydrostatic equilibrium, the codes 
must be able to model the cluster out to 1 Mpc while resolv- 
ing the 1 kpc core. This makes it a particularly challenging 
test. 

4.2 The Static Cluster 

4-2.1 Resolution of the cluster 

The two key requirements for success in this test are to be 
able to resolve the core and to have an accurate gravita- 
tional solver. Figure [TOl shows this case in point. It plots the 
density profile of the cluster at the start of the simulation 
and after 1 Gyr for steadily increasing levels of resolution. 
Although the King model does not have an analytical solu- 
tion, a one-dimensional numerical solution for the cluster's 
profile can be achieved from a simple numerical integration. 
This is shown in each plot as the solid line. In the top left im- 
age, the cluster is modelled on a static grid of size 100^ The 
initial match to the profile is good up to densities of 10~^, 
but the inner kiloparsec that contains the core is totally un- 
resolved. The effect is for the pressure to dominate over the 
gravity and the cluster expands, dropping the density in the 
core still further. After 1 Gyr, the cluster's profile is barely 
recognisable, despite the lack of external forces. Moving one 
plot to the right in Figure [TOl we see the results of adding a 
subgrid into the area that contains the highest density, i.e. 
the centre of the cluster. The density in the central region 
is now followed up to 0.15 Mqpc"''. The cluster is still not 
balanced and deviates away from the profile, but the shift is 
markedly less. As we continue to add in levels of refinement, 
we see the central density rise to meet the numerical expec- 
tation as the core becomes more resolved. The profile of the 
cluster changes progressively less until we reach five levels of 
refinement, when the change becomes almost undetectable 
except at the very centre. 



4-2.2 General results 

Figure [TT] shows the density, temperature and entropy pro- 
files for the cluster after 1 Gyr. For this test, Enzo was run 
with an initial grid of 100'' and eight additional levels of 
refinement, each reducing the cell size by a factor of two. 
These subgrids were placed anywhere where the cell mass 
was above a critical value. The cluster was set up in a box 
of size 3 Mpc with isolated gravitational boundary condi- 
tions. This gave a minimum cell size in the core of 0.12 kpc. 
Flash ran with the same boxsize and with periodic bound- 
ary conditions. It used a slightly larger initial grid of 128'' 
(since the gravity solver requires factors of two) and included 
seven additional levels of refinement, each of which reduced 
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Figure 10. Density profiles of tlie static King-model cluster run with Enzo (Zeus) with increasing resolution. The solid line shows the 
ID calculated numerical solution for the profile, with the red crosses and open blue boxes indicating the simulated profile at Gyr and 
1 Gyr, respectively. The number to the top right of each panel indicates the size of the root mesh and the number of refinements (i.e. the 
bottom right panel has eight levels of refinement each increasing in resolution by a factor of 2). 



the cell size by a factor of two. This gave a minimum cell 
size in the core of 0.18 kpc. Hydra and Gadget2 radially per- 
turbed a glass of 100,000 particles to the desired density 
profile within a periodic 3 Mpc box. 

As can be seen in the density profile in Figure [TTI with 
sufficient central resolution all the codes are successful at 
keeping the cluster in equilibrium and resolve the core well. 
All the codes show negligible deviations from the static King 
profile after 1 Gyr, matching the analytical profile over seven 
orders of magnitude, down to densities of 10"^". At the low 
density edge of the cluster, small deviations (~ 15%) from 
the profile are seen as the cluster edge diffuses into the back- 
ground. Enzo (PPM), Enzo (Zeus) and Flash produce the 
closest fit to the numerical solution at these densities, owing 
to the larger choice of initial grid which fixes the minimum 
resolution. Hydra and Gadgets show small deviations close 
to the cluster's tidal radius as they start to run out of par- 
ticles. 

Unlike the shock tests in section |3l all our codes in this 
test and in section refine based on increasing mass. Al- 



though grid based codes can make resolution choices from 
other parameters, mass in a grid cell (i.e. density) is the most 
used for simulations on scales above a few parsecs. However, 
in contrast to the particle based codes, a grid code's min- 
imum resolution is always known, since it is fixed by the 
specified minimum grid size. Depending on the simulation 
type, a fixed minimum resolution is either a huge advantage 
or a drain on computational time. In simulations where the 
denest structures are the primary focus, SPH's ability to 
sweep all the particles into these areas means that high res- 
olution can be gained quickly and efficiently. To gain the 
same efficiency from a grid code, a careful choice of mini- 
mum resolution grid and refinement criteria must be made. 
Both the test in this section and the following one are exam- 
ples of this problem type and were significantly more time 
consuming for grid based codes than for particle based ones. 
This is in contrast to the shock tests which ran faster. How- 
ever, if the focus of a simulation is the surrounding structure 
of a dense object, then the ability to specify a minimum res- 
olution is extremely useful and time efficient. In a particle 




Figure 11. Profiles from the static King-model cluster after 1 Gyr. Left to right shows density, temperature and entropy. The black 
line shows the expected profile. In this test, Enzo used an inital grid of 100'' and 8 levels of refinement (minimum cell size, Axmin, of 
0.12 kpc). Flash used an initial grid of 128^ with 7 levels of refinement (Ai'min = 0.18 kpc) and Gadget2 and Hydra used 100,000 particles. 



Figure 12. Density projections of the cluster over the course of 1 Gyr in which is moves once around the simulation box. Images taken 
at 0, 250, 500, 750, 1000 Myrs with projected density range [10* 10^®'^] Mq Mpc^^. Yellow and red shows higher density regions than 
green while black is very low density. (Images produced with Enzo (Zeus)). 



scheme this would normally be obtained by adding enough 
additional particles such that both the dense and surround- 
ing medium were resolved; a solution which may significantly 
add to the expense of the simulation. 

Figure [TT] demonstrates the two SPH codes achieve a 
similar result with 100,000 (46'^) particles to that obtained 
by the AMR codes with initial meshes of around 1 million 
cells and 8 or 9 levels of refinement. To achieve this with a 
static mesh would be impossible, as the mesh would need to 
be 25,600^ Not surprisingly then, given the complication 
of placing and interpolating between refinements the SPH 
codes are much faster. We note that in the SPH runs the 
distance from the centre to the 32nd particle is 0.4 kpc, 
which is very close to twice the minimum cell size of the 
AMR runs in the core. Again, we see that in order to obtain 
convergence between the AMR and SPH codes we require 
roughly one particle per cell in the region of interest. Un- 
fortunately for the AMR codes, for this particular problem 
which includes a wide density range the ability of SPH codes 
to naturally increase resolution with mass is very advanta- 
geous. 

Finally for this section, Figure[TT]also clearly shows that 
both of the Enzo implementations include noticeable ringing 
in the temperature profile which is reflected in the entropy 
and less obviously density profiles. This causes a deviation 



from the analytical solution of around ±20% in the inner 
2 kpc (10 cells). This ringing is not present in the analytic 
solution and doesn't appear in either Flash or the SPH codes 
which closely resemble both each other and the calculated 
solution. 

4.3 The Translating Cluster 

Using the stable clusters developed above we can test the 
Galilean invariance of the codes by giving them a velocity 
relative to the static simulation volume. This is a commonly 
encountered situation for cosmological simulations where 
large objects often move at many hundreds of kilometres per 
second relative to the background. At sufficiently high ve- 
locities it is well known and straightforward to demonstrate 
that mesh based codes are not Galilean invariant, w hereas 
particle based methods are (e.g. IWadslev et "al]|2008h . What 
we wish to investigate is the size of these effects for typical 
astrophysical objects moving at typical astrophysical veloc- 
ities. 

For this test the cluster was given a bulk velocity such 
that, in 1 Gyr, it moved around the simulation box once, 
returning to its original starting position. This is shown vi- 
sually in Figure[l2]which also demonstrates the box's bound- 
ary conditions which are periodic. Since there are no exter- 
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Figure 14. Profiles from the translating King-model cluster after 1 Gyr. Left to right shows density, temperature and entropy. The black 
line shows the expected profile. In this test, Enzo used an initial grid of 50^ and 9 levels of refinement (minimum cell size, Axmim of 
0.12 kpc), Flashused an initial grid of 128^ with 7 levels of refinement (Axniin = 0.18 kpc) and Gadget2 and Hydraused 100,000 particles. 



nal forces acting on the cluster, the end profile should be 
identical to the initial one. 

The time consuming nature of this test for grid codes 
also became an issue. Enzo (Zeus) and Enzo (PPM) lowered 
their initial grid to 50'^, but added an additional level of re- 
finement maintaining a minimum cell size of 0.12 kpc. Flash 
managed, after a long run, to maintain the same resolution 
as for the static test, with a minimum cell size of 0.18 kpc. 
The two SPH codes also had the same set-up as for the static 
case. 

Figure [13] shows the residue remaining if the initial con- 
figuration is subtracted from the final one. Another fun- 
damental difference between grid based and particle based 
techniques is visible here in the background medium. Un- 
like a particle code, grid codes cannot have zero density and 
energy cells. Therefore, a cool, low density gas surrounds 
the cluster in this test when run by both Enzo and Flash. 
In the case of Figure 1121 which only contains results from 
Enzo (Zeus), the minimum value shown is the background 
density, which is therefore displayed in black. From looking 
at the cluster core in Figure clearly neither of the Enzo 
simulations returns to the correct position, with a partic- 
ularly significant lag in the case of Enzo (PPM). In both 
these cases, the cluster as a whole correctly returns to the 
centre of the simulation box, but the central core does not, 
leaving it off-set from the centre of the cluster. This can 



be seen clearly in the last image of Figure 1121 After 1 Gyr, 
Enzo (PPM) has a core 176 kpc from the original start po- 
sition while Enzo (Zeus) finishes with its core 93 kpc short 
of the initial position. The outer regions of the cluster also 
show distortion due to the lower resolution of the initial grid. 
These problems don't appear to occur for either the SPH 
codes or Flash. Flash shows a small amount of distortion in 
the outer region of the cluster, although the central region 
stays largely intact with the core only 1.7 kpc from the centre 
of the box. Both SPH codes show very similar results with 
Hydra getting its core closest to the original start position 
with a 0.4 kpc off-set and Gadget2 with a 25 kpc off-set. In 
lower density regions, we see some distortion to the cluster's 
structure which becomes more marked at higher radii. This 
is due to the number of particles significantly decreasing as 
we move away from the core. 

Figure [14] is equivalent to Figure [11] except that now 
the cluster has moved once around the periodic box. All 
codes maintain the density profile of the cluster extremely 
well, with the only deviations appearing in the core. Here, 
Gadgets performs best, resolving the core in good agreement 
with the analytical prediction. All other codes struggle to re- 
solve the core after 1 Gyr, with Flash struggling the most and 
underestimating the core density by a factor of 10. Lower res- 
olution runs (minimum cell size of 0.23 kpc) using Enzo also 
showed a similar drop in core density in this test, whereas 
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a much smaller difference was noticed for the static clus- 
ter at the same resolution. This change in core density was 
seen, bu t not understood, in the Santa Barbara Comparison 
Project (|Frenk et alJI 19991 ) which observed the lower cluster 
densities in AMR codes. Without an analytical prediction 
for what the correct density should be, it was impossible to 
determine which of the recorded densities was more accu- 
rate. This test highlights two important effects controlling 
these core properties; resolution and advection. In the first 
instance, it is important to resolve the cluster core other- 
wise its size will be overestimated, leading to a lower central 
density. This is demonstrated clearest in Figure [TD] where in 
the case with four levels of refinement, the cluster remained 
intact but the central density still dropped. In the second 
case, the explicit movement of the fiuid through the grid (i.e. 
advection) in AMR can (but not necessarily does) cause a 
drop in the central density. This is an explanation for why 
the static and translating king models differ in the AMR 
case, even at the same resolution. We see from Figure [14] 
that Flash suffers most from advection issues whereas Enzo 
(PPM) shows minimal difference from the static model, fol- 
lowing the analytic profiles at least as well as the SPH codes. 
In both these areas, SPH's ability to follow the mass allows 
it to excel quickly and efficiently compared with AMR codes 
who have to be more careful when selecting their refinement 
criteria. 

While this test provides clues towards the deviation in 
cluster properties in the Santa Barbara Comparison Project, 
it should be noted that other effects can also be at work. 
In cosmological models, for inst ance, turbulence play s an 
important role as invest igated by I Wadslev et al.l (|2008l ) and 
[Mitchell et al.l (|in prep.l ). 

The temperature plot in the second panel of Figure 1111 
shows that Flash, Gadgets and Hydra all underestimate the 
core temperature whereas Enzo (PPM) and Enzo (Zeus) 
overestimate it slightly. Additionally, all the AMR code 
overestimate the cluster temperature at its edge with Enzo 
(PPM) having the poorest fit. The core's entropy is also 
overestimated in the case of the Flash and Enzo (Zeus); un- 
surprising since this depends on the inverse of the density. 

4-3.1 Gravitational tests: Gode-Specific Issues 

For both the static and translating cluster tests presented 
above, Enzo was run on a single processor. When using the 
parallised version of Enzo, the gravity solver suffers from ex- 
cess noise in the velocity motions which mix the gas, reduc- 
ing the entropy in the cluster centre. The result is a density 
drop in the cluster core and an expansion of its outer regions, 
unrelated to the issues discussed in section [4^ The result- 
ing density profiles using the parallised ve rsion of Enzo are 
shown in Figure [TSl In other situations fe.g. lTasker fc BrvanI 
L2OO8) the gravity solver has been tested successfully to give 
identical results in serial and parallel, but an improved, less 
noisy, gravity solver is currently being developed to solve 
this problem where is occurs. 

Enzo (PPM) also experienced problems dealing with 
the background medium. The sharp drop in conditions be- 
tween the cluster edge and low density background gas can 
cause Enzo (PPM) to artificially produce negative densities 
and energies. This was corrected for by setting a minimum 
density of O.lMoMpc"^, which negative densities were au- 




r (kpc) 

Figure 15. Density profile from the translating king model after 
1 Gyr run with Enzo in parallel. Run on multiple processors, the 
spurious velocity motions cause the gas to become over-mixed and 
the cluster density to drop in the core and expand at the cluster 
edge. The black line shows the expected profile. 

tomatically set to. Having to approximate densities multiple 
times during the run is a likely reason for the momentum 
loss that is clearly observed in Figure [141 

The default multigrid algorithm used in the publicly 
available version of Flash 2.5, proved too slow to be used 
practically. In order to overcome this we implemented a hy- 
brid FFTW based multigrid gravity solver, resulting in a 
vast reduction in the runtime of the simulations. The latest 
version of Flash 3. includes a more efficient multigrid algo- 
rithm and a nested FFTW gravity solver is currently under 
development for release in the near future (similar to the 
one used here). 

For Gadgets we encountered a problem in that although 
the radial profiles were very stable the cluster as a whole 
tended to drift around over time. This is due to the difficulty 
of very accurately determining the lowest order term in the 
gravitational force expansion for a treecode. Each individ- 
ual term includes a small error, with these errors largely but 
not exactly uncorrelated. Under these conditions the total 
momentum of the system is not guaranteed to be conserved 
exactly and a "random walk" occurs. As the configuration 
is designed to be entirely static the direction of the resid- 
ual force is highly correlated from one step to the next and 
despite recovering the correct value to one part in 10* this 
still leads to a net drift. It is possible to circumvent this, 
as shown above, by dramatically reducing the opening angle 
for the tree but this rapidly negates the advantage of using a 
tree in the first place. For more normal simulations this tiny 
zeroth order force error is of course negligible as the random 
motion of the particles disorders the direction of the drift er- 
ror as time progresses. P^M based gravity solvers like that 
employed by Hydra do not suffer from this problem as the 
zeroth order term in the Fourier transform is automatically 
zero. 

5 DISCUSSION AND CONCLUSION 

We have presented four well posed tests with known solu- 
tions which we have used to compare four well used astro- 
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physical codes; Enzo, Flash, Gadget2 and Hydra. We have 
examined each code's abihty to resolve strong shocks and 
incorporate an accurate gravity solver capable of resolving 
forces over many orders of magnitude in density. 

The tests presented here were specifically designed to be 
difficult and push the codes to their limit; appropriately so 
since many situations in astrophysics demand extreme phys- 
ical conditions. Despite these requirements, all the codes did 
pass our tests satisfactorily and what we highlight here is the 
strengths and weaknesses of the different techniques. 

Significant issues were Enzo (Zeus)'s failure to produce 
a spherical shock, a problem that leads to a lack of energy 
conservation and to poor recovery of the position of a Sedov 
blast. For extreme shocks the artificial viscosity implementa- 
tion of Gadget2 leads to spurious entropy driven bubbles and 
for very static configurations the centre-of-mass can slowly 
drift due to the difficulty of recovering the zeroth order term 
in the tree force. For a translating cluster, both implemen- 
tations of Enzo appear to lose momentum over time. 

All our tests lead to the conclusion that in order to 
achieve roughly the same resolution in grid and particle 
codes a good rule of thumb is to require one particle per 
cell in the region of interest. This is difficult to achieve for 
the SPIf codes over shock fronts or voids and likewise harder 
for AMR codes to match in the centre of collapsed objects. 
This is reflected strongly in our results where the AMR codes 
largely excel at the shock tests in section [3] and the SPH 
codes at the gravitationally bound cluster tests in sectional 
For a cosmological simulation these criteria would equate to 
having the minimum mesh cell size about the same size as 
the gravitational softening or half the SPH search length. 

As implied by the relation above, for strong shocks SPH 
codes require significantly more effort than AMR codes. 
This is because AMR codes have the ability to add extra 
resolution in regions of rapidly changing density, whereas 
SPH codes can effectively only refine with the density it- 
self, so they end up undertaking lots of unnecessary work 
far from the shock where the forces are small in regions of 
high but uniform density. However, for models such as the 
King sphere, where AMR codes traditionally refine using a 
density criterion anyway, SPH codes can achieving high res- 
olution far more efficiently in these dense structures because 
they do this very naturally by following individual particles. 
For AMR codes to achieve similar results, high resolution 
and a significant increase in computational time is required. 
In these cases AMR codes would still have an advantage if 
it was necessary to resolve the low density region, where the 
SPH codes simply run out of particles and have to smooth 
over very large volumes. 

Additional care must be taken for the parallel imple- 
mentations of Hydra and Enzo. Hydra's parallel version was 
not used in this paper and further testing would be needed 
to ensure that it performs correctly. The gravity solver in 
Enzo was demonstrated to produce different results in se- 
rial and parallel in the King model tests (although not in 
other situations) and this should be explicitly tested when 
running a simulation with this code. An updated gravity 
solver for Enzo which rectifies this problem is currently be- 
ing developed. While serial runs allow much parameter space 
exploration, rapid code development and inter-comparison, 
true state-of-the-art large simulations are confined to codes 



that can utilise high performance computing facilities and 
successfully run on a large number of processors. 

To a large extent the choice of simulation code largely 
comes down to the problem being attacked. For problems 
with large dynamic range or when the object is rapidly 
translating across the volume particle based methods are 
much less computationally expensive in order to achieve the 
same result. Conversely, for problems where the interesting 
physics is in regions of rapidly changing density rather than 
in high density regions themselves, AMR codes excel thanks 
to their more adaptive refinement criteria. 

The computational time required for each code to per- 
form these tests varied greatly. Ideally, elapsed time is 
largely irrelevant since the required resolution and problem 
size should be set by the appropriate physics being tackled 
rather than available computer resources. However, we in- 
clude Table[2]for completeness and to give the reader an indi- 
cation of times involved. In Table [2j the computational sys- 
tems referred to are; 'UF HPC is the University of Florida's 
high performance computing center which uses (typically) 
InfiniBand-connected 2.2 GHz AMD Opteron duel cores 
with 6 Gb RAM. 'UF Astro' is the University of Florida As- 
tronomy's minicluster with 2.66 GHz Intel quad-cores and 
4 Gb RAM. 'Columbia' is the Columbia University Astron- 
omy's cluster with gigabit-connected 2.2 GHz Opteron cores 
and 3 Gb RAM. 'Miranda' is at the University of Durham 
and consists of Myrinet-connected 2.6 GHz Opteron cores 
with 4 Gb RAM. 'Nottingham' is the Nottingham University 
HPC with gigabit-connected 2.2 GHz Opteron cores with 2 
Gb of RAM. 

While this paper attempts to cover the major fea- 
tures of astrophysical simulations, it does comprise only 
four tests. Further examples would extend this paper be- 
yond its scope (and readability), but the differences be- 
tween codes cannot be fully cataloged without further test- 
ing. This paper then, is designed as a starting point for 
a suite of tests to be developed from which codes can 
be quantitatively compared for the jobs they are intended 
for. To assist groups wishing to run these tests on their 
own code and to compare new updates, we are making 
these results and initial conditions available on the web at 
^http : //www. astro .ufl . edu/ codecomparison, 
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Table 2. Computational time to perform the presented tests, including the number of processors used and the type of system. Readers 
should note that comparing CPU times between simulations performed on disparate computers is fraught with danger. For true runtime 
performance, the codes would need to be bench marked on identical systems which was beyond the scope of this paper. We therefore 
include this table purely for guidance and warn against over interpretation, f The times quoted for the static and translating cluster 
runs in Flash are estimated total runtimes based on the speed up attained with the latest version of the FFTW gravity solver. The 
original runs were performed using an early version of the hybrid FFTW gravity solver which was still under development and as such 
was not optimised. The original runs took place on 24 processors and runtimes for the pre-development level FFTW solver were 217 
hrs 42 mins for the static cluster and 564 hrs 5mins for the translating cluster. This drastic improvements arises from a decrease in the 
memory requirements and improved communication patterns. A reduction in computational time in the Sedov blast test was additionally 
achieved by using a simplified Riemann solver that was not required to follow multiple fluids with different adiabatic indices (a feature 
the default configuration of Flash has implemented to model thermonuclear flashes). This reduced the memory overheads. 
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04:00 


1 


UF Astro (Intel Xeon) 


Flash 


00:24 


32 


Miranda (Opteron) 
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and Gadget2 simulations were carried out on the Notting- 
ham HPC facility. The Flash simulations were carried out 
on the Virgo Consortium computing facility in Durham. 
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